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D , Abstract 



We present results from a new Monte Carlo simulation for jet fragmentation in 
QCD and SUSY QCD for large primary energies ^/s up to 10 16 GeV. In the case of 
SUSY QCD the simulation takes into account not only gluons and quarks as cas- 
cading particles, but also their supersymmetric partners. A new model-independent 
hadronization scheme is developed, in which the hadronization functions are found 
from LEP data. An interesting feature of SUSY QCD is the prediction of a size- 
able flux of the lightest supersymmetric particles (LSPs), if R-parity is conserved. 
About 10% of the jet energy is transferred to LSPs which, owing to their harder 
spectra, constitute an important part of the spectra for large x = E/E- ]et . Spectra 
of protons and of secondary particles, photons and neutrinos, are also calculated. 
These results have implications for the decay of superheavy particles with masses 
up to the GUT scale, which have been suggested as a source of ultrahigh energy 
cosmic rays. 

Key words: Perturbative calculations in QCD, Supersymmetry, Cosmic rays 



1 Introduction 

QCD, being an essential part of the Standard Model, successfully describes accel- 
erator data for production of hadrons in e + e~ annihilation and in deep-inelastic 
scattering. There are two distinctive parts in these calculations: the perturbative 
QCD computation of the parton cascade in jets, and the parton hadronization, in 
which low-virtuality partons are converted non-perturbatively into hadrons. 

The QCD parton cascade is usually studied in Modified Leading Logarithmic Ap- 
proximation (MLLA), where large logarithms, ln(Q 2 ) and ln(l/:c), play a crucial role 
(here Q 2 is the maximum of the perpendicular momentum k±, and x = fci|/fc|^ ax ). 

This approximation is characterized by remarkable features. 

In the MLLA the QCD cascade has a probabilistic interpretation, provided by 
the absence of interference terms in the tree diagrams. The colour coherence effect 
is taken into account in the MLLA. It suppresses the emission of soft gluons and 
results in the Gaussian peak of the parton distribution in terms of £ = ln(l/x) 
(hump-backed plateau). 

The evolution of parton cascades in the MLLA (as well as in the LLA) is ad- 
equately described by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) 
equations |j]. 

The parton spectra can be obtained analytically and by Monte Carlo (MC) 
simulations. 

Examples for analytical solutions are the limiting spectrum [0] and the Gaussian 
spectrum Q, in which we include the distorted Gaussian spectrum [Q, ||. The 
limiting spectrum is the most accurate of them. Since we have a special interest in 
it, we shall shortly review below the basic assumptions under which this solution is 
obtained. 

The limiting spectrum gives the energy spectrum of partons -Du m (£,Y) for a 
given center-of-mass energy y/s of an e + e~-pair. Here D is the number of cascade 
partons, Y = ln(y / s/2A) and A is the dimensional QCD scale. 

The analytical expression for Aim(£> Y) is given in Refs. ||, ||. There are two 
fundamental parameters involved in the limiting spectrum solution: the scale A and 
the minimal virtuality Qq of partons, down to which the cascade develops pertur- 
batively; Qo can be viewed as the effective mass of the partons. Two assumptions 
are necessary for the validity of the limiting spectrum. The QCD coupling constant 
a s (k 2 _) evolves with k 2 ^ effectively as in the one-loop approximation with three fla- 
vors nj = 3 for all k 2 ^, 

/r2\ 127T 

as ^ )= (33-2n f )Hk 2 ± /Ai)- (1) 

As a matter of fact, A in Eq. (HI) is treated in the limiting spectrum solution as 
a free parameter to fit the e + e~-data. The best fit corresponds to A = 250 - 
270 MeV. For this range of A values, a s (Mz), given by Eq. (Ill) with n/ = 3 is in the 
interval 0.118-0.120, to be compared with the average experimental value a s (Mz) = 
0.1184±0.0031 |5J. Therefore, the phenomenological parameter A coincides well with 
A-QCD; which fits the experimental value of a s (Mz), in the one-loop approximation 
with nf = 3. 

The second assumption, necessary for the derivation of the limiting spectrum, 
is Qo = A. It gives a reasonable value of Qo, but the exact equality of these values 
has no theoretical justification. 



The limiting spectrum solution is valid only for small i< 1. In this region, which 
includes the maximal values of multiplicity (in the Gaussian peak), it describes very 
accurately experimental data at all available energies yfs (see, e.g., [||, |9|). The 
large x up to x = 1 give the dominant contribution to the total momentum of 
cascade partons. Therefore, the limiting spectrum solution does not guarantee that 
J xDn m (x, s)dx precisely equals 2, although it can be valid up to large x ||. 

Monte Carlo (MC) simulations of the QCD cascade give a more precise descrip- 
tion of the cascade evolution. They are valid for all x including x ~ 1. In contrast 
to the limiting spectrum, in MC simulations one can use a s (k^_) with the measured 
value of Aqcd, varying the number of flavours and two- loop corrections beeing 
taken into account. The assumption A = Q , specific to the limiting spectrum, is 
not needed. MC simulations are based on a probabilistic interpretation of the jet 
cascade. Parton branching is described by the Altarelli-Parisi functions, and the 
probability of parton evolution between two values of virtualities without branching 
is given by the Sudakov form factor. Finally, the coherent effect in the soft gluon 
emission (destructive interference) is conveniently taken into account by angular or- 
dering 9\ > 62 > 9^ . . . |jl0 |, where the indices number the generations (non-ordered 
processes are suppressed [ 1 1J) . The first MC simulation with angular ordering was 
performed in Ref. |l(J. At present there are several detailed MC simulations, e.g. 
[12, 13, [LJ, |15|| , which differ mainly in their description of hadronization. We shall 



now briefly discuss the problem of hadronization. 

The description of parton hadronization is based on the assumption of Local 
Parton-Hadron Duality (LPHD) fl6|| . This hypothesis implies that when Qq is 
small enough (of the order of A) there is a proportionality between the spectra of 
partons and hadrons, with relations between their momenta, which are local in the 
phase space. Such an interpretation can be justified by the idea of preconfinement 
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As far as spectra are concerned, LPHD implies a proportionality between the 
hadron and parton spectra. In Refs. [jj g], it is emphasized that, most reasonably, 
this proportionality holds not on a "one parton - one hadron" basis, but for the 
number of particles averaged over a finite interval A£ ~ 1. 

The LPHD hypothesis for limiting spectrum straightforwardly results || in 

Dhad(x, Vs) = K h (Q )D pa _ n (x, y/s, Q ), (2) 

where the constant K^ is universal, in the sense that it does not depend on ^/s. 

Equation (0) completes our description of the limiting spectrum, expressing 
the hadron spectrum through the spectrum of partons. The constant Kh, which 
connects the two spectra, is found from a comparison with experimental data as 
Kh ~ 1.3 for A = Qq ~ 270 MeV ||, and it does not change with energy unless 
some new physics (e.g. supersymmetry) appears. 

In MC simulations, the parameter Qq is in principle a free parameter found 
by fitting experimental data. For HERWIG ||] and PYTHIA [||], for example, 
Qo ~ 1 GeV. Several detailed hadronization models are used in simulations, e.g. 
the independent fragmentation model, the Lund string model fl~8||, and the cluster 



fragmentation model [19|. Usually, these models use many free parameters and 
require to keep track of the four-momentum evolution of all partons. 

The calculations described above are valid up to y/s ~ 1-10 TeV. At higher en- 
ergies the production of supersymmetric particles is expected to change the results. 



One might be interested in much higher energies, being inspired by the production 
of superheavy particles up to the GUT scale in the Universe. Such particles can be 
produced by topological defects and by many processes at the post-inflationary stage 
of the Universe. Recently superheavy particles with masses Mx ~ 10 12 -10 14 GeV 
attracted much attention as a source of the observed Ultra High Energy Cosmic 
Rays (UHECR) with energies 10 19 -10 20 eV (for recent reviews, se e plf ). 

The limiting spectrum for SUSY QCD was calculated in Ref. J22J for very high 
energies y/s, corresponding to masses of superheavy particles Mx ~ 10 12 -10 14 GeV. 
The super symmetric partons (squarks and gluinos or jointly spartons) participate 
in the cascade until the virtualities t of the particles drop below the mass scale of 
SUSY particles, t ~ -^fusY- Then a SUSY particle decays, producing in the end 
the Lightest Supersymmetric Particle (LSP), for which the lightest neutralino is 
usually considered. The role of supersymmetric partners is two-fold: they double 
the number of parton types in the cascade, and they change the evolution of a s (k'j_). 
Even at small t <C -^susy> the casc& de remembers the number of flavours at large 
t because, for example, each squark leaves after its decay a quark, which continues 
QCD cascading. At large t and small x <C 1 gluons and gluinos dominate and their 
"children" constitute the dominant part of the cascade at small t. Therefore, the 
dominant contribution to the limiting spectrum is given by gluons and gluinos. 

The SUSY QCD limiting spectrum solution has two drawbacks with respect to 
ordinary QCD. First, the number of flavours that determine the evolution of the 
coupling constant according to Eq. (||) has to be fixed to one value of nj for the 
whole range of k 2 ^. Second, the limiting spectrum for ordinary QCD is normalized 
by experimental data, which are absent in the case of SUSY QCD. Normalization 
due to the conservation of momentum f xDn m (x, s)dx = 2 is unreliable, since the 
limiting spectrum is not valid for large x, which give the main contribution to the 
integral (see discussion in [p2]] ). 

During the last few years, the production and decays of supersymmetric particles 
have been included in most MC simulations focusing on LHC studies. Although the 
LHC will operate above the expected threshold of SUSY particle production, its 
energy is not large enough for these particles to participate in the QCD cascade. 
Therefore, all currently available MC simulations consider only on-shell decays of 
spartons and neglect possible branchings of gluinos and squarksR. Another obstacle 
against the use of standard MC simulations at extremely large energies around 
y^ ~ 10 12 — 10 14 GeV is that the necessary numerical precision and the required 
amount of memory space and computing time become a challenge for present-day 
computers. 

We have therefore developed a new MC simulation, which includes as cascading 
particles not only gluons and quarks but also gluinos and squarks. We consider 
cascades that are initiated by the decay to two jets of superheavy particles with 
mass Mx ~ 10 12 -10 14 GeV, or by e + e~ annihilation at s = M x . SUSY partons, 
squarks and gluinos, are produced in the fragmentation of ordinary partons and 
vice versa. All squarks and gluinos are assumed to have equal masses, for which 
we use Msusy = 200 GeV and M S tjsy = 1000 GeV. When the virtually of the 
cascading particles drops below M| USY , spartons decay to LSPs (neutralinos), which 
freely escape. The perturbative development of the cascade continues with ordinary 
partons until their virtualities reach Q 2 ,, for which we use Qq = 0.625 GeV 2 to 



-'The future C++ version of HERWIG will include branchings of spartons 



fit the data at small energies y/s. We use a new hadronization procedure. It is 
based on a model- independent, phenomenological approach, in which hadronization 
functions for large s (or Mx) are calculated from hadron spectra observed at small s 
(Mjf). This method can be used for any type of hadrons as well as for photons and 
neutrinos, if their spectra are known with good enough accuracy at small energy. 

Following Ref. f20f| , we shall use the following notation: 
A is the dimensional QCD scale, 
Y = ln(V5/2A), 

t = p 2 is the virtuality of cascade partons, 
Q 2 = tmax is the virtuality of the primary parton, 

Qq is the minimum virtuality of the perturbative evolution of the QCD cascade, 
z = E'/E, where E and E' are the energies of ingoing and outgoing partons at 
fragmentation, 

£ = 1 — cos 8, where 6 is the angle between two outgoing partons, 
i=CE\ 

fej_, feii are the transverse and parallel momenta transferred, respectively, 
x = fe||/fep x , 
£ = ln(l/x). 



2 MC Simulation of the Perturbative Phase 
of SUSY QCD Cascades 

The perturbative part of our simulation is very similar to those of MCs for ordinary 
QCD cascades, except for including spartons and the condition for their exit from 
the cascade. We consider a superheavy X particle with mass Mx which decays into 
two jets with energy E- ]et = Mx/2. We assume that the primary partons produced 
in the X particle decay have the maximum virtuality Q 2 = m^/4 and that the 
X particle has equal branching ratios to all partons. As to the first assumption, 
in reality, there is a distribution of partons with different t, but the Sudakov form 
factors suppress small t values. The second assumption is made because of the 
unspecified interactions of the X particles. 

Our simulation closely follows the angular ordered parton cascade algorithm 
developed in Refs. |l(], [19| ]. It is convenient to use in this algorithm the variable 
i = C,E 2 , where E is the energy of the incoming parton, Q ~ 1 — cos#, and 6 is 
the angle between the two emitted partons. A primary parton with energy E- ]e t 
(= rax/2) and angular variable £o < 1 initiates a cascade, which proceeds until 
the ordinary partons reach the minimal virtuality t = 4Qq- Here the perturbative 
evolution of the cascade terminates. 

In each branching of an incoming parton i with t', we generate with the veto 
algorithm p4j a new t and z according to the probability distribution 



dtdz 
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dV l (t,z) = J2 T ^-» s z 2 {\-zft P^ ik (£)-±>. (3) 
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Here, the sum includes all possible branching channels jk, z 2 (l — z) 2 t is the parton 
transverse momentum, and Aj is the product of the individual Sudakov-like form 



factors Ai^jk 

Ai^ jk (i) = exp 
with 



At ■ t 



t n-+jk(t ) 



(4) 



— a s [z 2 {l-zft\Pi^ jk {z). (5) 

-min 

The unregularized Altarelli-Parisi splitting functions Pi^jk(z) of SUSY QCD |26| 
are given in Table |]. 

The angular ordering (j,(k < Ci f° r the branching i — > j'/c, which takes into 
account colour coherence, is equivalent to t,- < z 2 £j and i& < (1 — z) 2 ii. These 
conditions result in 

■2min — V »min/f > ^max — •!• y '-min/'-- (,Oj 

For the evolution of the running coupling a s as a function of gluon virtuality t at 
small momentum transfer t < £susy> we use the standard two- loop dependence with 
variable nf and thresholds, and normalize a s as a s {Mz) = 0.119, which corresponds 

(5) 
to A-^j^ = 222 MeV. At large momentum transfer t > ^SUSY we use minimal SUSY- 

SU(5) coupling constant evolution J25|], normalizing the coupling constant at yt = 
M GU T = 1 • 10 16 GeV, as a s (M,2, UT ) w 1/25.8. Explicitly we use 

n m = »(^gut) m 

sU l + 6 s /(4vr)ln(i/M2 UT )a(M2 UT )' [) 

where b s = 9 — nj is a constant, that governs the evolution of the coupling constant 
with t. At £ > tsusY) n f = 6 and 6 S = 3. The above assumption means that 
we introduce, instead of many thresholds corresponding to SUSY particles with 
different masses, a single threshold at t = tsiJSY- This is a reasonable thing to 
do in view of the large uncertainties in our knowledge of mass spectrum of SUSY 
particles. Equatio n ([/] ) approximates accurately enough the evolution of a s (t) as 
calculated in Ref. [p7| 7 when isusY ~ 2 • 10 5 GeV 2 . Starting from this value, a s (t) 
evolves in the regime of Eq. (|7|). Note that isiJSY does not necessarily coincide with 
the scale Msusy> the universal mass of squarks and gluino, for which we use as two 
representative values Mstjsy = 200 GeV and Msusy = 1 TeV. In particular, the 
low value of £susy used here is compatible with much larger Mstjsy > as emphasized 
in Ref. ||. 

Finally, we have to specify the value of the cut-off i m j n for the cascade evolution. 
We do not distinguish between different quark flavours and use t m { n = 0.625 GeV 2 
for all branchings in which only ordinary particles are produced and i m \ n = M^gy , 
where MgusY is the typical mass scale of the spartons, for branchings in which 
SUSY particles are produced, respectively. 

Let us now describe a step i — ► jk in our simulation. For an incoming parton 
i with t' we generate first a new cascade variable t, according to the probability 
distribution given by the ratio Aj(i')/Aj(t). Then we select the branching channel 
jk using as weight fi—tjk{t) and generate z according to the probability distribution 
a s [z 2 {l - zft\ Pi^k(z). 

The last ingredient in the perturbative part of our simulation is the exit of 
supersymmetric particles from the cascade. We assume that the neutralino x 1S the 



LSP and that R-parity is conserved. Reaching i m - m = M^gy, squarks and gluinos 
decay as q — ► q + \ an d 9 ~^ Q + Q + X> thus producing UHE LSPs. 

In fact, we are running in this work two Monte Carlo similations: with the 
ordinary QCD and with SUSY QCD. 

In the former case supersymmetric partons are not included, and for perturbative 
calculations we assume the SM particle content with a s (t) evolution in two-loop 
approximation with proper thresholds. We fix Qq = 0.625 GeV 2 . We need these 
calculations mostly to test our method. 

The assumptions of SUSY QCD Monte Carlo are described above. At i min < 
Mg USY cascade develops according to ordinary QCD scheme. 



3 Hadronization 

The Monte Carlo simulation described in the last section is completely determined 
by perturbative physics. How the spectrum of coloured quarks and gluons Di(x, i/s) 
is transformed into the spectrum of hadrons -Dhad(x> V^s) is still an open problem. 
Monte Carlo simulations have to use some hadronization model (see Introduction), 
which describes the non-perturbative evolution of the cascade for t < 4£ m j n . Two 
hadronization models, the cluster fragmentation model |1£] used in HERWIG [12] 
and the Lund string model [18] used in PYTHIA [13], require the knowledge of 
the four- momenta of all the partons. Thus these models need detailed time and 
memory-consuming computations. 

We suggest here a phenomenological, model-independent hadronization scheme 
based on the knowledge of the hadron spectra at energies y^ smaller than the 
energy of interest. This method is valid for any hadron type and can be applied 
to the secondary particles, such as photons and neutrinos, as well. The application 
of this method is somewhat restricted (e.g. it cannot give the angular distribution 
of particles in a jet or correlations), but its use is very efficient for the decay of 
superheavy particles, where multiplicity, and hence the number of partons to follow 
in a simulation is very large. 

Our hadronization scheme depends on only one theoretical assumption, which 
is reliable and testable. Namely, we assume that the unknown non-perturbative 
physics can be factorized into hadronisation functions fi(z) that do not depend on 

D h (x, V~s) = £ f 1 - D % {x/z, ^)ft(z) , (8) 

J x Z 

where the index h runs through different types of hadrons, e.g. tt ,^, N etc. 

The functions fi{z) give the probability that a parton i with energy E is con- 
verted into a hadron h with energy zE. It is implicitly assumed in Eq. @ that the 
perturbative cut-off Qq is fixed, and / 4 - is determined for this value of Qo, although 
in principle for every Qq and Di(x,yfs,Qo) one can find f^(z,\/s,Qo) to fit the 
observed hadron spectra. 

Equation (g) with energy independent hadronization functions follows from basic 
principles and is confirmed (see below) at energies of e + e~ colliders. It has the form 
of a Volterra integral equation of the first kind though in contrast to the standard 
case, the RHS contains not one but two unknown functions fi for every h. In 
principle, the two functions f g (x) and f q (x) can be uniquely determined if D^ is 



known as an analytic function without errors for two different values of y/s. In 
practice, D^x) is known only as a discrete set of experimental data and Eq. (pi) 
represents an ill-posed inversion problem^ |pGJ]. Instead of solving Eq. (||) by an 
inversion method, we prefer to find physically motivated trial functions for /j to fit 
the experimental data at y/s = 91.2 GeV. 

In terms of the more convenient variable £ = ln(l/x), Eq. (J8|) has the form 

Dh& v^) = E / d ^ A(e - e', >/i)/i(0 , (9) 

i=q,9 ° 

where the index h in the hadronization functions is suppressed. 

In the limiting spectrum, when Qq = A, the hadronization functions fy are 
proportional to delta functions. Inspired by this analytical solution we choose for 
fi Gaussian functions, 

fi(0 = a i exv(- ^-^ 2 \ . ( 10 ) 

With this hadronization function the approximate proportionality holds between 
spectra of partons and hadrons as LPHD demands. The position of the peak in 
the hadronization function determines the shift between the maxima of parton and 
hadron spectra. While for gluons the hadronization function / 9 (£) should vanish 
for £ — ► 0, because gluons have to split their energy to a qq pair, for quarks / 9 (£) 
can be finite at £ = 0. 

The hadronization functions we obtained for Qq = 0.625 GeV 2 from a fit to LEP 
data at y/s = 91.2 GeV are shown in Fig. [j]. 

Our hadronization scheme has been tested by two methods: for relatively small 
energies, y/s = 58 GeV and 133 GeV, we confronted our calculations with LEP 
data, and for very large y/s (or Mx) we compared the calculated spectrum with the 
limiting spectrum, using a special case when it is correct (see below). In both cases 
ordinary QCD Monte Carlo was used. 

Figures 2-4 display a comparison between the charged hadron spectrum from 
our MC simulation for ordinary QCD and experimental data [28] at y/s = 58, 91.2 
and 133 GeV, respectively. 

Let us now discuss whether the hadronization functions /i(£) found from the fit 
to data at y/s = 91.2 GeV can really be used at M x = 10 12 ~10 16 GeV. 

First of all we note, that a test can be given by LPHD, which demands approx- 
imate proportionality between parton and hadron spectra. It implies that the £' 
values, that give the dominant contribution to the integral in Eq. (|9]), are about the 
same at y/s = 91.2 GeV and at large Mx- Numerical tests show that this is indeed 
the case for both the quark and the gluon contribution. 

As a critical test of our hadronization scheme, we compared the limiting spec- 
trum with the results of our simulation for a special, well-controlled case of ordi- 
nary QCD with the number of quark flavours nj = 3, and with a s (/c 2 _) given by 



2 Volterra integral equations of the first kind can be solved normally by linearization, even if the LHS 
are data. However, the lower integration limit in (||) does not represent a sharp cut-off because the kernels 
Di(x) vanish for x — > 1. Therefore, Eq. (ph behaves effectively like a Fredholm equation, and these are 
known to be extremely ill-conditioned. 



Eq. (|TJ) . For the limiting spectrum in this case we can use the normalization constant 
Kh w 1.3 obtained by fitting experimental data Q. Since we do not introduce new 
high-energy physics, the limiting spectrum is valid for any initial energy Mx (see 
Introduction). In Fig. 0, we show the ratio of these two (charged) hadron spectra. 
The agreement between the two spectra is excellent, except for the small £ < 6 
region where it is known that the limiting spectrum is not valid. The disagreement 
reaches 50% at f ra 2.1 (x « 0.12). 

In conclusion, we think that our hadronization recipe is a valid alternative to the 
extrapolation of the Lund string or the cluster fragmentation model to extremely 
large Mx- 



4 Results: Spectra of Hadrons and Secondary 
Particles 

Using the algorithm for the perturbative evolution of the SUSY QCD cascade as 
described in Section 2 and our hadronization scheme from Section 3, we can now 
compute the fragmentation spectra of hadrons. As numerical values for Mx, we 
choose in the graphs given as examples three values interesting for UHECR physics, 
M x = 10 12 , 10 13 , 10 14 GeV, as well as M x = 10 6 and 10 16 GeV as lowest and highest 
scale of interest. Similarly, we use Msusy = 200 GeV and Msusy = 1000 GeV as 
two representative values for the SUSY mass scale. 

In Figs. |8|and|9], the hadron spectra dN^^/d^ from SUSY QCD MC simulations 
are displayed as function of £ for Msusy = 200 GeV and Msusy = 1000 GeV, respec- 
tively; in both figures the spectra were calculated for Mx = 10 12 , 10 13 , 10 14 GeV. 
For the GUT scale Mx = 10 16 GeV, the hadron spectra dN^^/d^ are shown in 



Fig. [TO] and for the low scale Mx = 10 GeV in Fig. 11. The hadron spectra de- 
pend only weakly on Msusy > with increasing differences for larger values of Mx- 
Both effects are easy to understand: when spartons disappear from the cascade at 
t ~ Mg USY due to on-shell decays, each of them leaves there an ordinary parton 
with similar virtuality. Therefore, the cascade proceeds as if nothing had happened, 
except that some energy is lost through the emission of neutralinos and leptons, 
which is not large (~ 10%). Second, the importance of spartons for the cascade 
decreases for smaller values of Mx, thereby reducing also the dependence of the 
hadron spectra on Msusy for smaller Mx ■ 

The weak dependence of the spectra on Msusy justifies our choice of an universal 
value for the masses of supersymmetric particles. Indeed, if we assume now that 
supersymmetric particles have different masses in the range 200-1000 GeV, the 
resulting hadron spectra will differ less than in Figs. |lJJ and [ll]. 

The signature of supersymmetry in decays of superheavy X particles is the pro- 
duction of LSPs, which we assume as stable neutralinos. They are generated in the 
cascade mostly when the virtuality of the spartons approaches M^ SY . The calcu- 
lated neutralino spectra are shown in Figs, p^-15 for the same parameters as the 



hadron spectra in Figs. Ef-jlj]. Like the hadron spectra, they have the characteristic 
Gaussian form, however with a shifted position of their maxima due to their larger 
cut-off Msusy in the shower development. The energy fraction taken away by the 
neutralinos is typically 10% for values of Mx interesting for UHECR physics, with 



a minimum of 5% for Mx = 10 6 GeV and MgusY = 1 TeV and a maximum of 12% 
for M x = 10 16 GeV and M SUSY = 200 GeV. 

We have only derived a common hadronization function for all hadrons and, 
consequently, we cannot calculate directly, e.g. pion or nucleon spectra through 
Eq. (|8|) . Since the fraction of energy e^ going into different meson and baryon species 
is determined by the non-perturbative process of hadronization, these fractions as 
the hadronization functions themselves do not depend on s. Thus, we can use the 
value from Z decay, ejv ~ 0.05 and e^ rs 0.95. Then 

diVnuci diV had dN n _ c dJV had 

dx dx dx ^ dx 

Using the hadron spectra obtained in the last Section, it is simple to calcu- 
late analytically the spectra of secondary particles, photons and neutrinos. The 
normalized photon spectrum from a decay of one X particle at rest is given by 

diV T 2 fids dAWi 



dx sHlf^f' <12> 

The total neutrino spectrum, given by the sum from decays of pions and muons, 
can be presented in the following form, 

dN u 2 (dN u .dJV, s, dN »e, ,\ c,o\ 

where for pion decay 

JAT l-l A„. A AT. . 

(14) 





dN Vli f 1 dy diV had 
(vr -► pup) =R 
dx Jrx y dy 


and for muon decay 




dx^ 


,. .. c) - P f ld y f y/r d y' dN »* diV had 
M Jx y J y y' dy dy' 


with 

dN„ e _ 2 

dy 


, dN u „ 5 9 4 , 
-6y 2 + 4y 3 , -^ = - - 3y 2 + -y 3 


and r = (m^/m^) 2 , R = 


:l/(l-r). 



(15) 



(16) 



The resulting nucleon, photon and neutrino spectra x dNi/dx are shown as func- 



tions of x together with the spectra of neutralinos in Fig. |17| for Msusy = 200 GeV 
and Mx = 10 12 GeV and Mx = 10 14 GeV, respectively. We have multiplied the 
spectra by x 3 in order to facilitate the comparison of our spectra with the en- 
ergy spectra of observed UHECR. At x > 0.7, the spectra have some uncertainties 
because of the unknown branching ratios of the X particle into (s)partons and fluc- 
tuations due to the small number of produced particles. The excess of nucleons over 
secondary particles, neutrinos and gamma, at largest x is a result of kinematical 
effect and the steep spectra of pions and nucleons in the end of the spectrum. 



5 Discussion 

In this Section we compare for large Mx the results of our MC for the two cases, 
SUSY QCD and ordinary-QCD, with other computations, and most notably with 
limiting spectrum calculations. The latter case of ordinary QCD is formally a special 
case of our MC simulation for SUSY QCD in the limit Msusy> ^susy — ► oo, i.e. a s 
is given by two- loop approximation with variable rif, and the probability to produce 
a sparton is zero. 

The validity of our method has been proved by the tests described in Section 
3. If no new physics beyond the three light quark flavours is introduced, the k±- 
dependence of a s is given by Eq. (|l]) with nj = 3 and the limiting spectrum with 
Kh = 1.3 is valid for arbitrary high energies. We can calculate the hadron spectrum 
in our ordinary QCD MC (hadronization procedure included), introducing there 
the same assumptions about n/ and a s . The excellent agreement is illustrated by 
Fig. 5. The disagreement seen at large x is natural, because the limiting spectrum 
is not valid there. 

It is instructive to compare our MC for ordinary QCD with variable rif and the 
exact behaviour of a s (k±) with the limiting spectrum with fixed number of flavours 
rif = 3 and rif = 6. It is clear that, in neither case, a s (k±) from Eq. (||) describes 
correctly a s in the whole interval of kj_, and the MC spectrum should be between 
these two solutions. Figures. 6 and 7 show that this is indeed the case. The accuracy 
of each limiting spectrum compared with the MC spectra is better than 30-50%. 

In Ref. |3(]], HERWIG was used to obtain fragmentation spectra in case of ordi- 
nary QCD. The maximal mass Mx possible to simulate was Mx = 10 11 GeV and 
even for this not very large value of Mx the computations required several months. 
The spectra were displayed only for large x > 0.01, beyond the Gaussian peak. One 
of the conclusions of this work was that at x > 0.2 the proton yield is higher than 
the photon and neutrino yield. However, it was later realized that this result is 
caused by the tendency of HERWIG to overproduce protons at large x (Ref. |3l| , 
see also [J33J] ). 

Let us come over to our SUSY QCD MC and compare the simulated spectra 
with the SUSY limiting spectrum |22|| . The spectra disagree both in the position 
of the Gaussian peak and in its height. To clarify which assumptions of the SUSY 
QCD limiting spectrum are responsible for this disagreement, we re-run the SUSY 
QCD MC simulation with a set of assumptions as similar as possible to those used 
in the derivation of the SUSY QCD limiting spectrum. We found that the main 
reason for the disagreement is the universal dependence of a s (t), taken as aj 1 (t) = 
(6 s /4vr)ln(t/A 2 ), with b s = 3 for SUSY, together with A = Q = 250 MeV. It 
differs from a s with a variable number of flavours, which is used in SUSY QCD 
MC, by a factor 1.4-3 in the whole k\ range, with largest disagreement at small 
k±. Changing the evolution of a s (k±), an agreement can be reached between the 
MC and the SUSY QCD limiting spectrum: we run the SUSY QCD MC including 
only gluons and (massless) gluinos with fixed b s = 3 and with frozen a s (t) for 
i < 0.9 GeV 2 , which is a reasonable physical assumption. The comparison with the 



SUSY QCD limiting spectrum for partons is shown in Fig. 16. The two spectra 
agree indeed quite well. 

An interesting alternative approach to computing the fragmentation spectra pro- 
duced by decays of superheavy particles was suggested in a recent work |H|] . In this 
method, the event generator SPYTHIA [32] was used to simulate fragmentation 



spectra of partons and spartons into protons, photons, and neutrinos at the scale 
Mx = 10 4 GeV. Then the DGLAP equations were used to evolve the fragmentation 
functions up to the scale 10 12 -10 13 GeV. 

It is premature to compare our results, since in [HTJ preliminary results are pre- 
sented, but the spectra, as displayed in [^] and |33|], do not agree well with ours. In 
particular, the Gaussian peak is broader than in our calculations, Fig. |l7]. Compar- 
ing these spectra one should be aware of the differences in methods and assumptions. 
For example, we treat spartons as cascading particles, while in SPYTHIA spartons 
are taken as on-shell particles, which decay but do not cascade. On the other hand, 
the SPYTHIA spectrum is used only as the input, and the evolution to higher en- 
ergies includes cascading. This difference will be eliminated with the C++ version 
of HERWIG H|, which wil1 be available soon. 



6 Summary 



We have developed a new MC simulation for jet fragmentation in ordinary QCD and 
SUSY QCD, which is valid for initial energies up to the GUT scale. The simulation 
includes a perturbative part, operating at virtualities higher than the infrared cut-off 
Qq = 0.625 GeV, and a hadronization part. 

The perturbative part for SUSY QCD includes squarks and gluinos as cascade 
particles with a universal mass Mstjsy- The evolution of a s (t) takes into account 
the correct number of active flavours and spartons at a given t. The influence 
of the scale Mstjsy on the hadron spectrum is rather weak for the studied range 
300 < Mstjsy < 1000 GeV. It implies that if supersymmetric particles have different 
masses in the range 200 -1000 GeV, the resulting difference in the hadron spectra 
remains small (see Section IV and Figs. |l0| and |TT| ). 

The hadronization scheme is model-independent and based on the well justified 
and tested assumption that the hadronization function /^(z), see Eq. (||), does not 
depend on y/s. Thus the hadronization function could be calculated from LEP data. 
Our scheme was tested at yfs = 58 and 133 GeV against experimental data and for 
very large values of Mx by a comparison with the limiting spectrum solution for 
ordinary QCD. For the aim of this comparison, we calculated hadron spectra using 
the MC for ordinary QCD in the case when the limiting spectrum is known to be 
correct: rif = 3 and a s (t) given by Eq. (|l|). The excellent agreement between both 
spectra is illustrated by Fig. ||. 

The spectra of nucleons and secondary particles, photons, neutrinos, as well 



as neutralinos, have been calculated and presented in Fig. |17|. These spectra can 
be used for calculations of fluxes of ultra high energy cosmic rays, produced by 
superheavy dark matter and by topological defects. 
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Table 1: Splitting functions Pi-+jh(z), where z is the energy fraction of the particle j. 
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Figure 1: Hadronization functions for quarks f g (£) (solid line) and for gluons f g (£) 
(broken line) obtained by fitting Gaussians to experimental data at \fs = 91.2 GeV. 




Figure 2: Comparison of the spectrum of charged hadrons diV c h/d£ from ordinary QCD 
Monte Carlo simulation (solid line) with the experimental data (shown with errorbars) at 
yfi = 58 GeV. 
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Figure 3: Comparison of the spectrum of charged hadrons diV c h/d£ from ordinary QCD 
Monte Carlo simulation (solid line) with the experimental data (shown with errorbars) at 
Vs = 91.2 GeV. 
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Figure 4: Comparison of the spectrum of charged hadrons diV c h/d£ from ordinary QCD 
Monte Carlo simulation (solid line) with the experimental data (shown with errorbars) at 
y/s = 133 GeV. 




Figure 5: The ratio R = -Diim(£)/-^Mc(£) of the limiting spectrum and of the hadron 
spectrum from the simulation for Mx = 10 12 GeV (solid line), 10 13 GeV (broken line) 
and 10 14 GeV (dashed line). All for ordinary QCD with nf = 3. 
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Figure 6: Comparison of the limiting spectrum of QCD for rif — 3 and rif = 6 with the 
ordinary QCD spectrum from the Monte Carlo simulation: R = Z?ii m (QCD, n/ = 3) /Due 
(solid line) and R = D hm (QCD,n f = 6)/D MC (broken line). Both for M x = 10 12 GeV. 
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Figure 7: Comparison of the limiting spectrum of QCD for n/ = 3 and rif = 6 with the 
ordinary QCD spectrum from the Monte Carlo simulation: R = Aim(QCD,n/ = 3)/Dmc 
(solid line) and R = Aim (QCD, n/ = 6)/D MC (broken line). Both for M x = 10 14 GeV. 







Figure 8: Hadron spectra d/V^/d^ from SUSY QCD Monte Carlo simulation for Mx 
10 12 GeV (bottom), M x = 10 13 GeV (middle) and M x = 10 14 GeV (top), all for M SUS y 
200 GeV. 
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Figure 9: Hadron spectra diV/j/d£ from SUSY QCD Monte Carlo simulation for Mx 
10 12 GeV (bottom), M x = 10 13 GeV (middle) and M x = 10 14 GeV (top), all for M SU sy 
1 TeV. 
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Figure 10: Hadron spectra cW/j/d£ for SUSY QCD Monte Carlo simulation for Msusy 
200 GeV (broken line) and M SU sy = 1 TeV (solid line), both for M x = 10 16 GeV. 
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Figure 11: Hadron spectra dNh/dC, from SUSY QCD Monte Carlo simulation for 
Msusy = 200 GeV (broken line) and M SU sy = 1 TeV (solid line), both for M x = 10 6 GeV. 




Figure 12: Neutralino spectra from SUSY QCD Monte Carlo simulation for M x 
10 12 GeV (bottom), M x = 10 13 GeV (middle) and M x = 10 14 GeV (top), all for M SUS y 
200 GeV. 
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Figure 13: Neutralino spectra from SUSY QCD Monte Carlo simulation for Mx 
10 12 GeV (bottom), M x = 10 13 GeV (middle) and M x = 10 14 GeV (top), all for M SU sy 
1 TeV. 
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Figure 14: Neutralino spectra from SUSY QCD Monte Carlo simulation for Msusy 
200 GeV (top) and M SU sy = 1 TeV (bottom), both for M x = 10 16 GeV. 
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Figure 15: Neutralino spectra from SUSY QCD Monte Carlo simulation for Msusy 
200 GeV (top) and M SU sy = 1 TeV (bottom), both for M x = 10 6 GeV. 







Figure 16: Parton spectrum from SUSY QCD Monte Carlo simulation (boxes) and of 
the SUSY QCD Limiting Spectrum (solid line) for Mx = 10 12 GeV. Both for gluons 
and gluinos only. The coupling constant in the Monte Carlo simulation is frozen at 
t < 0.9 GeV 2 , and b s = 3 is fixed in both cases. 
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Figure 17: Neutrino, gamma and nucleon fragmentation spectra from SUSY QCD Monte 
Carlo simulations for Mx = 10 12 GeV (solid lines) and 10 14 GeV (dotted lines), all for 
Msusy = 200 GeV. 



